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Foreword 


This report represents work completed in the 1997-1998 time frame while working as 
a NASA contractor for SAIC at NASA Langley Research Center. It was originally 
submitted as a paper to Smart Materials and Structures in Feb., 1998 (Log Number: 
00395) and conditionally accepted in July, 1998 pending submission of a minor revision. 
This was withheld by the author because it was simultaneously submitted as part of a 
patent disclosure, LAR-15827, in April, 1998. 

Much has changed since the time this report was written. This model has since been 
duplicated and improved upon over the past 9 years since it was written. It is being 
reproduced here for historical purposes and to make it more referenceable for those who 
may wish to reference it. 
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Abstract 

An analytic approach is developed to predict the shape and displacement with voltage 
in the quasi-static limit of LaRC Thunder Actuators. The problem is treated with classical 
lamination theory and Von Karman non-linear analysis. In the case of classical 
lamination theory exact analytic solutions are found. It is shown that classical lamination 
theory is insufficient to describe the physical situation for large actuators but is sufficient 
for very small actuators. Numerical results are presented for the non-linear analysis and 
compared with experimental measurements. Snap-through behavior, bifurcation, and 
stability are presented and discussed. 

Introduction 

Thunder actuators are devices constructed from an isotropic laminate of aluminum, 
LaRC SI adhesive, piezo-electric PZT, and a metal backing of either steel or brass. The 
LaRC adhesive is a 1 mil thick solid thermoplastic and the aluminum layer is indented to 
provide electrical conductivity after the actuator is cured at 250 deg C. After the curing 
process they form a bow shape after reaching room temperature. If simply supported or 
arranged in a clamshell arrangement they produce a linear motion when an external 
voltage is applied. 

A simple method for controlling the manufacturing process and predicting the 
actuation amplitude has long been sought for Thunder actuators. A first attempt to do this 
by an analytic method was completed by BBN [9]. However, this theory treated the 
problem in 1 dimension with no adhesive layer (a three layer model) and it wasn't clear if 
this was sufficient to describe the behavior. In particular, there are many variables to be 
considered. To make assumptions regarding the shape of the device can be misleading. 
Quite different shapes are possible depending on the geometry and the material 
properties. The other behavior missing from the BBN model is the snap-through 
behavior. These devices have two natural bending modes - one in the x and one in the y. 
In reality, this behavior is a strong function of size. Below some critical size the shape is 
a perfect dome with equal curvatures but above this critical size the curvatures become 
unequal. 
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One approach to this problem is to use classical lamination theory. This is a very old 
subject based on the theory of composite laminates. With laminates, temperature effects 
are very important. As a result a very rich theory has been developed to minimize these 
temperature effects. Although the theory is sufficient to find an optimal lay-up to 
minimize temperature effects, it is quite deficient in predicting the actual shape in all 
cases. The reason for this is because it is a linear approximation to a problem that is 
essentially non-linear. As a result, a non-linear correction to the theory is needed to 
describe large actuators. 

The type of correction used here is one that is valid for large displacement and small 
strain. This correction is usually attributed to Von Karman [2] in plate theory. Hyer [3] 
has applied this to his treatment of orthotropic laminates with good success. He uses this 
in combination with a simple Raleigh-Ritz energy method to obtain solutions exhibiting 
the snap-through behavior. To apply his theory here, however, it must be reformulated for 
isotropic laminates, which is a completely separate (and unrelated) case. In particular, the 
qualitative behavior of orthotropic laminates and isotropic laminates are quite different. 
Orthotropic laminates produce saddle shapes that approach cylinders when the laminate 
is large. Isotropic laminates produce domes that approach cylinders when the laminate is 
large. They are both subsets of the anisotropic case but not subsets of each other. This 
behavior is well demonstrated with this theory. 

Solution by Classical Lamination Theory 
Derivation of equations 

Classical lamination theory [1] assumes strains are small, and amplitudes are shallow. 
In addition, it assumes both shearing and extensional strains perpendicular to the surface 
are zero. These assumptions follow the Kirchhoff-Love hypothesis for shells. 

Under deformation the thin plate is described by three variables - u, v and w, where u 
is the displacement in x direction, v in the y and w is the height of the plate in the z 
direction. The variables u, v and w are themselves functions of x, y and z. It is assumed u 
and v follow the relation 


u = u n - z 


v = v n 


d_% 
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3 w 0 
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where the 0 subscript denotes the quantity is measured about the middle surface. The 
strains are given by [2] 
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The results of Equation 1 and Equation 2 imply 
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where k x . Ky, and K X y are the x, y and twist curvatures. The middle surface strains and 
curvatures are given by. 
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The relationship between stress and strain is given by 
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where Q is the reduced stiffness and Q is the transformed reduced stiffness. For isotropic 
media. 
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where Y is Young's modulus and v is Poisson's ratio. To find the force and moments we 
first find the stresses for the k'th layer by. 
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Now we use the relationship 
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where h is the thickness. By Equation 9 and 8 we find, 
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where A, B, and D are the extensional, coupling and bending stiffness matrices and are 
given by, 
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D «= 3 l.feWS -*-')• 

where the z's are measured with respect to z=0 (at a point half way through the 
thickness). Since each layer is isotropic the symmetry of the stiffness matrices will be the 
same as the reduced stiffness matrix. Moreover, since in this case the layers are not 
symmetric about z=0 the coupling stiffness will be non-zero so that, 
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We now have everything we need to find the thermal and actuation loading so we 
may solve the problem. 

Thermal Loading 

The thermal loading takes two forms. The first is extensional and the second is 
bending. The problem is greatly simplified since we have only isotropic layers with equal 
thermal strains in each direction. The thermal loading may be calculated by, 


where, 
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Actuation Loading 

By analogy with thermal loading, 
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where t is the thickness of the PZT, V is the applied voltage, d 3 i is the actuation constant 
in a plane perpendicular to the poling axis and, 


N A =d 3 fa+Q^V, 
M A =d 3l (Q l + Ql(z 3 +$V. 


(16) 


Here it is assumed we only have actuation on the third layer. 

Thermal and Actuation Displacement 

To find the resulting strains and curvatures in terms of the actuation loading we first 
add both the thermal and actuation sources as 
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Now we use the results of Equation 10 so that 
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Inversion of the above gives 
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By using the results of 12 and inverting the 6x6 matrix we find, 

‘ ” (A 1 + A,)(fl 1 +D 2 )-(S 1+ S,) 2 ’ 


(A, + A 2 ) M - ( B ! + /V 

(a 1 + a 2 )(a+a)-(A+5,) 2 ’ 


(19) 


( 20 ) 


Y .xy ~ ^ XV ~ 0 • 

So we have thus found the curvature in terms of known quantities. The shape is that 
of a dome. The simplicity of the relationship comes about because the layers are 
isotropic. The solution is valid for small actuators such as the round ones used for pumps 
and speakers. 


Solution by the Von Karman Non-Linear Approximation 

The Von Karman approximation assumes large displacement and small strain. This 
suggests including second order terms in w in addition to the linear terms for u, v and w. 
Expanding Equation 2 to include these terms give, 
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This combined with the results of Equation 1 give, 


e = 


0 w 0 1 

0 x 2 


0 W 0 

0 X 


- z 


0 " W 0 

~dx r ’ 


£ y = 


0 v 0 , 1 (dm 


0 y + 2 l 0 y 


- z- 


d-w 0 
0 y 2 


( 22 ) 


_ 0 Up , 0 Vp , 0 W 0 dw 0 " 0 2 w 0 

^ 03 ; dx dx dy " 0 x dy 

The total energy takes the form, 

U = / / / (dU 0 -dW T -dW^, (23) 

where dU 1 is the stored elastic energy volume density and dWj and dWA are the 
temperature and actuation contributions. The elastic energy portion of it for the k'th layer 

is. 
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For the actuation portion we have 


(dW^ = (o A x e x + a$ e) 3 , 
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(26) 


Now that we have the general form for the energy integral we may try a Raleigh-Ritz 
type solution. The method involves making a guess for the solution. This guess should 
look qualitatively like the expected solution yet be flexible enough to allow for 
adjustment. This adjustment comes in the form of changing the guess based on a finite 
number of parameters that the guess contains. The parameters are adjusted in such a way 
as to minimize the total energy of the system. In this situation it would be advantageous 
to choose a guess that is close enough to the classical lamination result that in the linear 
limit (small scaling) the solution approaches the classical lamination result. Such a guess 
takes the form 


w° = j (a x 2 + b y 2 ). (27) 

Classical lamination theory assumes no shearing strain between layers for thermal 
expansion. If we make the same assumptions here we find u*-* and v^ must take the form, 

u° = c x- ^ a 2 x } - ^ ab x y% 

(28) 

v° = d y-X a 2 f - \ a b x 2 y. 
o ‘ 4 

The ansatz represented by Equation 27 and 28 has been used successfully by Hyer [3] 
in his treatment of orthotropic laminates. It is sufficiently flexible enough to be used with 
isotropic laminates as well as we will soon see. Equations 22, 27 and 28 give, 

s x =c-X a bf-az, 

E y =d-^ab r-b z, (29) 

Y.vy = 0 - 

All that is needed now is to use Equation 23 with Equations 24, 25 and 26 and 29 to 
find the total energy. The solution is obtained by minimizing the total energy with respect 
to a, b, c and d. The result of the minimization is 
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The A's, B's and D's have the same meaning as before and N' is the same as well. The 
lengths L x and Ly are the length and width of the actuator. It is a simple matter to show 
that the solution in the linear limit (small scaling - L x =Ly=0) that the solution to these 
equations matches the classical lamination result if we take a=-K x , b=-Ky, c=s x and d=Sy. 
The non-linear equations are coupled and third order so it is unlikely that an exact 
analytic solution exists. It may be possible to find close approximate solutions to them, 
however. In any event they may be solved numerically with no difficulty. 
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Unlike the classical lamination result the non-linear analysis predicts unequal 
curvatures that depend on the magnitude of the scaling and the aspect ratio. It also 
predicts multiple possible solutions. The nature of these solutions are dome-like solutions 
which approach cylinders in the limit of large scaling. 

As a comparison for a particular case we compute a simple example. For this case 1 
mil aluminum, 1 mil LaRC Si, 5 mil PZT and 2 mil stainless was chosen. 

The non-linear analysis predicts a total of 5 possible solutions. Out of the five, 3 are 
non-physical so we are left with two real solutions. The snap-through behavior is clearly 
exhibited in this case. As a comparison this is compared with classical lamination theory. 
The behavior is quite different in this case. It should be noted that the non-linear theory 
may admit very different solutions depending on the geometry, aspect ratio, etc. 


-o.o: 

-o.o 

-o.c 

-O.i 


c 



l 


Figure 1. Non-linear Solution #1 
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1 

Figure 2. Non-linear Solution #2 


Stability 

The stability of each solution is determined by whether or not the solutions are 
maxima or minima. If a particular solution represents a minimum then it is stable. If it is 
a maximum, it is unstable. The total energy of the system may be viewed as a function of 
several variables. As such, once an extremum is found, the test whether it represents 
minima is determined by whether the determinant of F is positive definite where 

Det[F] > 0, 

(32) 

F - d 2 U 
' j d Xi d x j 

Fortunately d U/d xi is represented by the left hand side of Equation 30 (to within a 
multiplicative factor) so all we have to do are the xj differentiations to find the matrix 
elements. 

Bifurcation 

As stated earlier, The solution may take on different shapes depending on the scaling. 
To demonstrate this, numerical results were completed for a square piece having the same 
material properties as before except we vary the size from x=0 to x=2 in. For a small 
piece the shape is a dome as predicted by classical lamination theory. At some critical 
size the shape turns suddenly into a dome with unequal curvatures. In the limit of large 
scaling it becomes a cylinder. Figure 3 is a plot of both the x and y curvatures. In this 
case the critical size is at approximately 1 inch square. This is quite different from the 
results of Hyer[3], who found saddle solutions. However, this is easily explained by the 
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fact that his work was with orthotropic laminates as opposed to isotropic laminates. In the 
case of orthotropic laminates even classical lamination theory predicts a saddle. 



Figure 3. Bifurcation plot of x and y curvature as a function 
of size for a square actuator. Stable solutions only. 


Preliminary calculations seem to indicate the currently manufactured devices are 
close enough to the bifurcation point for these effects to be important. 

Experimental Results 

As an example we compare the model with known measured data, (assumed) material 
properties, and thicknesses. A table of these samples is shown below. 




Sample 1 (2"xl"). 




Material 

Thickness 

Y (GPa) 

CTE( 10-6/deg C) 

Layer 1 

Aluminum 

.001" 

73 

22 

Layer 2 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 3 

PZT 5a 

.01" 

63 

1 

Layer 4 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 5 

stainless steel 

.005" 

200 

16 




Sample 2 (2"xl"). 




Material 

Thickness 

Y (GPa) 

CTE( 10-6/deg C) 

Layer 1 

Aluminum 

.001" 

73 

22 

Layer 2 

LaRC Si 

.001 "(irregular) 

4 

46 
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Layer 3 

PZT 5a 

.01" 

63 

1 

Layer 4 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 5 

stainless steel 

.003" 

200 

16 



Sample 3 (2"xl"). 




Material 

Thickness 

Y (GPa) 

CTE(10- 6 /deg C) 

Layer 1 

Aluminum 

.001" 

73 

22 

Layer 2 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 3 

PZT 5a 

.01" 

63 

1 

Layer 4 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 5 

stainless steel 

.001" 

200 

16 



Sample 4 (3"xl.75"). 




Material 

Thickness 

Y (GPa) 

CTE( 10-6/deg C) 

Layer 1 

Aluminum 

.001" 

73 

22 

Layer 2 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 3 

PZT 5a 

.0078" 

63 

1 

Layer 4 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 5 

1020 steel 

.005" 

200 

16 



Sample 5 (3"xl.75"). 




Material 

Thickness 

Y (GPa) 

CTE( 10-6/deg C) 

Layer 1 

Aluminum 

.001" 

73 

22 

Layer 2 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 3 

PZT 5a 

.0078" 

63 

1 

Layer 4 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 5 

1020 steel 

.003" 

200 

16 




Sample 6 (3"xl.75"). 




Material 

Thickness 

Y (GPa) 

CTE( 10-6/deg C) 

Layer 1 

Aluminum 

.001" 

73 

22 

Layer 2 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 3 

PZT 5a 

.0078" 

63 

1 

Layer 4 

LaRC Si 

.001 "(irregular) 

4 

46 

Layer 5 

1020 steel 

.001" 

200 

16 


We compare the measured dome height with the predicted values. The first 3 samples 
were taken from the BBN report and the last 3 were measured at LaRC. The setting 
temperature of LaRC SI was assumed to be 250 deg C for the calculation. 
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Table 1. Measured vs. Predicted dome heights 


Sample 

1 

2 

3 

4 

5 

6 

Measured 

.114" 

.126" 

.193" 

.198" 

.281" 

.466" 

Predicted 

.116" 

.103" 

.029" 

.29" 

.282" 

.113 

Sol. Type 

cylinder 

cylinder 

dome 

cylinder 

cylinder 

cylinder 

BBN pred. 

.098" 

.095" 

.047" 

* 

* 

* 


The results of Table [1] show that the model is fairly accurate except when the base 
layer has a small thickness. This would suggest the adhesive layer irregularities have an 
effect on the solution. In the case of the smaller BBN 1 mil sample a transition from 
cylinder to dome was predicted. This would suggest the solution is close to the 
bifurcation point. Since any small error near the bifurcation point produces huge 
differences in the curvatures, this could be another explanation for the discrepancy for a 
small base layer. A third explanation could be the second order non-linear correction may 
not be enough to accurately describe the real physical shape in all cases. It is interesting 
the BBN model predicts the same trend. The reasons for this are not clear since they do 
not consider the adhesive layer at all. 


Conclusions 

More work needs to be done to accurately describe the problem to take into account 
irregularities in the adhesive layers. These irregularities are currently part of the 
manufacturing process and are currently needed to provide electrical conductivity for the 
device. At least in the case of the top layer it may be possible to treat the 
aluminum/adhesive combination as a single layer with special material properties of it's 
own . In the case of the steel/adhesive boundary, detailed measurements need to be taken 
to see if there is some trend in the irregularities. If there is it may be possible to model it. 
In any event, the model appears to be very accurate as it is for a base layer of 3 mils and 
above. 

This model shows some very interesting behavior and is a giant leap forward in 
demonstrating the snap-through behavior. It also gives insight into how the devices 
should be manufactured. It shows that in some cases perfect domes are possible and that 
one should construct the devices as far from the bifurcation point as possible so that the 
device is in the cylinder range. 

Numerical work needs to be done to show the voltage effects in Equation 16. In 
addition, external loading effects need to be taken into account. This is a very 
straightforward problem using the energy approach and would result in a slight variation 
of the Equation 30. Several variations of this are possible. The first would be to find the 
change in amplitude with voltage as a function of load. Another more interesting case 
would be to constrain the amplitude and allow the load to vary. This would result in an 
equation for blocked force as a function of constrained distance. 


16 



Acknowledgements 

I would like to thank Professor D. S. Cairns of Montana State University for his many 
very helpful comments and suggestions. In addition, I would like to thank Bill Hunter, 
Stephanie Wise, Tony Jalink, Bob Copeland, Robert Fox and Richard Hellbaum (all of 
LaRC) for helping provide measurements and other information. 

References 

1. Mechanics of Composite Materials 

R. M. Jones 
McGraw-Hill 1975 

2. Anisotropic Plates 

S. G. Lekhnitskii 

Gordon and Breach Publishers, 1968 

3. M. W. Hyer 

Calculations of the Room-Temperature Shapes of Unsymmtric Laminates 
J. Composite Materials 
Vol. 15, July 1981 

4. A. Hamamoto and M. W. Hyer 

Non-Linear Temperature-Curvature Relationships for 

Unsymmetric Graphite -Epoxy Laminates 

Int. J.Solid Structures 

Vol 23, No. 7, pp 919-935, 1987 

5. M. W. Hyer and P. C. Bhavani 

Suppression of Anticlastic Curvature in Isotropic and Composite Plates 

Int. J. Solid Structures 

Vol. 20, No. 6, pp 553-570, 1987 

6. R. Yang and M. A. Bhatti 

Non-linear Static and Dynamic Analysis of Plates 
Journal of Engineering Mechanics 
Vol. Ill, No. 2, Feb., 1985 

7. D. A. Pecknold and J. Ghaboussi 
Snap-Through and Bifurcation in a Simple Structure 
Journal of Engineering Mechanics 

Vol. Ill, No. 2, Feb., 1985 


17 



8. J. F. Campbell and D. S. Cairns 

A Multi-Purpose Sensor for Composite 
Laminates Based on a Piezo-Electric Film 
J. Composite Materials 
Vol. 26, No. 3 /1992 

9. Alan R.D. Curtis 

A Model of the THUNDER Actuator 

BBN Tech Memo #1190 

NASA Contract NAS 1 20101, Dec. 31, 1997 

10. Morgan Matroc Limited 

Piezoelectric Ceramic Products (product manual) 


18 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 
Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 
shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 


5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 


5f. WORK UNIT NUMBER 

810031.07.03 


11. SPONSOR/MONITOR'S REPORT 
NUMBER(S) 

NASA/TM-2007-2 14872 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category 37 

Availability: NASA CASI (301) 621-0390 

13. SUPPLEMENTARY NOTES 

An electronic version can be found at http://ntrs.nasa.gov 

14. ABSTRACT 

An analytic approach is developed to predict the shape and displacement with voltage in the quasi-static limit of LaRC 
Thunder Actuators. The problem is treated with classical lamination theory and Von Karman non-linear analysis. In the case of 
classical lamination theory exact analytic solutions are found. It is shown that classical lamination theory is insufficient to 
describe the physical situation for large actuators but is sufficient for very small actuators. Numerical results are presented for 
the non-linear analysis and compared with experimental measurements. Snap-through behavior, bifurcation, and stability are 
presented and discussed. 


15. SUBJECT TERMS 

PZT; THUNDER; Von Karman; Actuator; Bifurcation; LaRC-Si; Nonlinear; Piezo-electric; Rayleigh-Ritz; Snap-through 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 

L-19355 

10. SPONSOR/MONITOR'S ACRONYM(S) 

NASA 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

NASA Langley Research Center 
Hampton, VA 23681-2199 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


5d. PROJECT NUMBER 
5e. TASK NUMBER 


3. DATES COVERED (From - To) 
5a. CONTRACT NUMBER 


1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 

01-05 - 2007 Technical Memorandum 

4. TITLE AND SUBTITLE 

Quasi-Static Analysis of LaRC THUNDER Actuators 


6. AUTHOR(S) 

Campbell, Joel F. 


16. SECURITY CLASSIFICATION OF: 

17. LIMITATION OF 
ABSTRACT 

18. NUMBER 
OF 

19a. NAME OF RESPONSIBLE PERSON 

a. REPORT 

b. ABSTRACT 

c. THIS PAGE 

PAGES 

STI Help Desk (email: help@sti.nasa.gov) 






19b. TELEPHONE NUMBER (Include area code) 

u 

u 

u 

uu 

24 

(301) 621-0390 


Standard Form 298 (Rev. 8-98) 

Prescribed by ANSI Std. Z39.18 



































